library(tidyverse)  # ggplot(), %>%, mutate(), and friends 
library(broom)
library(MatchIt)  # Match things
library(Rcpp)
library(MASS)
library(modelsummary)
library(IRdisplay)
library(haven)
library(dplyr)
library(jtools)
library(car)
library(margins)
library(marginaleffects)
library(cobalt)
library(ggplot2)
library(scales)

# set directory


# IMPORT MAP
map <- read_csv('C:/./dataset76_20.csv')

regions <- read_csv('C:/./states_region.csv')
regions <- dplyr::rename(regions, State=state_x)


# replace string for being identical acrsso map and nibrs regions
map$State[map$State == 'Rhodes Island'] <- 'Rhode Island'
map$State[map$State == 'District of Columbia'] <- 'District Of Columbia'
# adding missing states
regions <- regions %>% add_row(...1= 888, State = "Alaska", country_region = "west", country_division="Pacific")
regions <- regions %>% add_row(...1= 999, State = "California", country_region = "west", country_division="Pacific")
regions <- regions %>% add_row(...1= 45878, State = "Florida", country_region = "south", country_division="South Atlantic")
regions <- regions %>% add_row(...1= 97666, State = "New Jersey", country_region = "north east", country_division="Middle Atlantic")


# add region and division to map dataset
map<- merge(map, regions, by="State")
library(janitor)

#clean labels: can be done by simply
map <- clean_names(map)

# change to character exposure variable
map$vic_race_black <- as.character(map$vic_race_black)

# final cleaning
map <- subset(map, map$agentype != 4)
map <- subset(map, map$state != "PAPSP8")

##################### IMPORT NIBRS
nibrs <- readRDS('C:/./nibrs_prepared_dataset_91_20.rds')

nibrs <- clean_names(nibrs)
nibrs$race_of_victim_black <- as.character(nibrs$race_of_victim_black)
nibrs$incident_date <- as.POSIXct(nibrs$incident_date) # otherwise margins won't work

# final cleaning
nibrs$agency_indicator[nibrs$agency_indicator == '5'] <- "other agencies"












#################################################################
####################### MAIN MODELS #############################
#################################################################



############## NIBRS (all with agency-level info and states)
m.out3ni <- matchit(race_of_victim_black ~ total_offender_segments+
                      total_victims + 
                      sex_of_victim+
                      agecat+
                      decade+
                      agency_indicator+
                      monthly_agency_overlap +
                      state_x,
                    data=nibrs,
                    method = "exact",
                    estimand="ATT")

bal.tab(m.out3ni)

new.names <- c(state_x_Idaho = "State: Idaho",
               state_x_Utah = "State: Utah",
               agency_indicator_city = "Agency: City",
               state_x_Montana = "State: Montana",
               agency_indicator_state_police = "Agency: State Police",
               monthly_agency_overlap = "Monthly Agency Homicide Overlap",
               sex_of_victim_male = "Victim's Sex: Male",
               sex_of_victim_female = "Victim's Sex: Female",
               "state_x_New Hampshire" = "State: New Hampshire",
               state_x_Oregon ="State: Oregon",
               state_x_Colorado="State: Colorado",
               "state_x_South Dakota"="State: South Dakota",
               "state_x_West Virginia"="State: West Virginia",
               "state_x_North Dakota"="State: North Dakota",
               state_x_Washington="State: Washington",
               state_x_Vermont="State: Vermont",
               state_x_Arizona="State: Arizona",
               "state_x_New Mexico"="State: New Mexico",
               state_x_Hawaii="State: Hawaii",
               state_x_Oklahoma ="State: Oklahoma",
               "agecat_21-25"="Victim's Age: 21-25",
               "agecat_71-75"="Victim's Age: 71-75",
               "agecat_76-80"="Victim's Age: 76-80",
               "agecat_81-85"="Victim's Age: 81-85",
               state_x_Massachusetts="State: Massachusetts",
               "agecat_61-65"="Victim's Age: 61-65",
               state_x_Michigan="State: Michigan",
               "agecat_16-20"="Victim's Age: 16-20",
               "agecat_66-70"="Victim's Age: 66-70",
               state_x_Iowa="State: Iowa",
               state_x_Nebraska="State: Nebraska",
               state_x_Ohio="State: Ohio",
               "agecat_56-60"="Victim's Age: 56-60",
               agency_indicator_other_agencies="Agency: Other Agencies",
               state_x_Kansas="State: Kansas",
               "state_x_South Carolina"="State: South Carolina",
               "agecat_51-55"="Victim's Age: 51-55",
               "agecat_26-30"="Victim's Age: 26-30",
               state_x_Maine="State: Maine",
               state_x_Kentucky="State: Kentucky",
               "agecat_0-5"="Victim's Age: 0-5",
               total_victims="N of victims",
               "agecat_46-50"="Victim's Age: 46-50",
               "agecat_86-90"="Victim's Age: 86-90",
               agency_indicator_tribal="Agency: Tribal",
               state_x_Virginia="State: Virginia",
               state_x_Texas="State: Texas",
               "agecat_41-45"="Victim's Age: 41-45",
               total_offender_segments="N of offenders",
               state_x_Tennessee="State: Tennessee",
               "agecat_91-95"="Victim's Age: 91-95",
               state_x_Wisconsin="State: Wisconsin",
               state_x_Missouri="State: Missouri",
               state_x_Georgia="State: Georgia",
               "state_x_North Carolina"="State: North Carolina",
               "state_x_Delaware"="State: Delaware",
               "agecat_6-10"="Victim's Age: 6-10",
               state_x_Alabama="State: Alabama",
               "agecat_31-35"="Victim's Age: 31-35",
               decade_90s="Decade: 1990s",
               "state_x_Rhode Island"="State: Rhode Island",
               state_x_Louisiana="State: Louisiana",
               decade_10s="Decade: 2010s",
               state_x_Nevada="State: Nevada",
               "agecat_36-40"="Victim's Age: 36-40",
               state_x_Arkansas="State: Arkansas",
               state_x_Mississippi="State: Mississippi",
               state_x_Maryland="State: Maryland",
               state_x_Connecticut="State: Connecticut ",
               "agecat_96-99"="Victim's Age: 96-99",
               state_x_Illinois="State: Illinois",
               decade_00s="Decade: 2000s",
               state_x_Wyoming="State: Wyoming",
               "state_x_New York"="State: New York",
               state_x_Minnesota="State: Minnesota",
               state_x_Indiana="State: Indiana",
               state_x_Pennsylvania="State: Pennsylvania",
               "agency_indicator_covered-by another agency"="Agency: covered by oth. agency",
               "state_x_District Of Columbia"="State: DC ",
               "agency_indicator_university or college"="Agency: University/college",
               "agecat_11-15"="Victim's Age: 11-15"
)

p<- love.plot(m.out3ni, stats = c("mean.diffs"),
          thresholds = c(m = .1, v = 2), abs = TRUE, 
          binary = "std",
          var.order = "unadjusted",
          var.names=new.names,
          colors=c("sienna3", "deepskyblue3"),
          sample.names=c("Unmatched", "Matched"),
          title=NULL)

p + scale_x_continuous(labels = function(x) {
  ifelse(x < 1, gsub("^0", "", format(x, nsmall = 2)), format(x, nsmall = 2))
})

matched_dataset3ni <- match.data(m.out3ni)



###### modeling 

#adjusted
model_matched3ni <- glm(solved ~
                          race_of_victim_black + total_offender_segments+total_victims + 
                          sex_of_victim+
                          agecat+
                          decade+
                          weapon+
                          circumstance+
                          agency_indicator+
                          factor(monthly_agency_overlap) +
                          state_x+
                          population_group+
                          location_type,
                        family=quasibinomial(), data = matched_dataset3ni, weights=weights)

tidy(model_matched3ni)





# get odds ratios and 95% cis
summ(model_matched3ni, exp = TRUE, robust=TRUE, digits=3)



# GENERAL ATT
ame_3ni<-comparisons(model_matched3ni, variables = "race_of_victim_black", vcov="HC3",
                     wts="weights",
                     newdata = subset(matched_dataset3ni, race_of_victim_black == 1)) |>
  summary()
ame_3ni$dataset <- 'NIBRS (1991-2020)'
ame_3ni$model <- 'Adjusted/Matched'

# simulation-based 
sim_model_matched3ni <- clarify::sim(model_matched3ni, vcov = "HC3")
est_sim_model_matched3ni <- clarify::sim_ame(sim_model_matched3ni, var="race_of_victim_black", subset=race_of_victim_black==1, contrast="diff", verbose=TRUE)
plot_est_sim_model_matched3ni <- plot(est_sim_model_matched3ni)+theme_classic()+geom_density(fill="blue")
xlab('Estimation')+facet_wrap(vars(.data$est),labeller=as_labeller(c('Ciao', 'Ciao2', 'Ciao3')))



summary(est_sim_model_matched3ni)



# DECADE ATT
ame_3ni_dec <- comparisons(model_matched3ni, variables = "race_of_victim_black", 
                           by="decade", vcov="HC3", wts="weights",
                           newdata = subset(matched_dataset3ni, race_of_victim_black == 1)) |>
  summary()

ame_3ni_dec$dataset <- 'NIBRS (1991-2020)'
ame_3ni_dec$model <- 'Adjusted/Matched'

ame_3ni_dec_pw <- comparisons(model_matched3ni, variables = "race_of_victim_black", 
                              by="decade", vcov="HC3", wts="weights",
                              newdata = subset(matched_dataset3ni, race_of_victim_black == 1),
                              hypothesis="pairwise") |>
  summary()

ame_3ni_dec_pw$dataset <- 'NIBRS (1991-2020)'
ame_3ni_dec_pw$model <- 'Adjusted/Matched'


# SEX ATT
ame_3ni_sex <- comparisons(model_matched3ni, variables="race_of_victim_black", 
                           by="sex_of_victim", vcov="HC3", wts="weights",
                           newdata = subset(matched_dataset3ni, race_of_victim_black == 1))|>
  summary()

ame_3ni_sex$dataset <- 'NIBRS (1991-2020)'
ame_3ni_sex$model <- 'Adjusted/Matched'

ame_3ni_sex_pw <- comparisons(model_matched3ni, variables="race_of_victim_black", 
                              by="sex_of_victim", vcov="HC3", wts="weights",
                              newdata = subset(matched_dataset3ni, race_of_victim_black == 1),
                              hypothesis="pairwise")|>
  summary()

ame_3ni_sex_pw$dataset <- 'NIBRS (1991-2020)'
ame_3ni_sex_pw$model <- 'Adjusted/Matched'





######### unmatched
model_unmatched3ni <- glm(solved ~
                            race_of_victim_black + total_offender_segments+total_victim_segments + 
                            sex_of_victim+
                            agecat+
                            decade+
                            weapon+
                            circumstance+
                            agency_indicator+
                            factor(monthly_agency_overlap) +
                            state_x+
                            population_group+
                            location_type,
                          family=binomial, data = nibrs)

summ(model_unmatched3ni, exp = TRUE, robust=TRUE, digits=3)
# GENERAL ATT
ame_3ni_unmatched<-comparisons(model_unmatched3ni, variables = "race_of_victim_black",
                               newdata = subset(nibrs, race_of_victim_black == 1)) |>
  summary()
ame_3ni_unmatched$dataset <- 'NIBRS (1991-2020)'
ame_3ni_unmatched$model <- 'Adjusted/Unmatched'

# simulation-based 
sim_model_unmatched3ni <- clarify::sim(model_unmatched3ni)
est_sim_model_unmatched3ni <- clarify::sim_ame(sim_model_unmatched3ni, var="race_of_victim_black", subset=race_of_victim_black==1, contrast="diff", verbose=TRUE)
plot_est_sim_model_unmatched3ni <- plot(est_sim_model_matched3ni, color='red')+ylab("Densità")




summary(est_sim_model_unmatched3ni)

# DECADE ATT
ame_3ni_unmatched_dec <- comparisons(model_unmatched3ni, variables = "race_of_victim_black", 
                                     by="decade", newdata = subset(nibrs, race_of_victim_black == 1)) |>
  summary()

ame_3ni_unmatched_dec$dataset <- 'NIBRS (1991-2020)'
ame_3ni_unmatched_dec$model <- 'Adjusted/Unmatched'

ame_3ni_unmatched_dec_wp <- comparisons(model_unmatched3ni, variables = "race_of_victim_black",
                                        by="decade", 
                                        hypothesis="pairwise", newdata = subset(nibrs, race_of_victim_black == 1)) |>
  summary()

ame_3ni_unmatched_dec_wp$dataset <- 'NIBRS (1991-2020)'
ame_3ni_unmatched_dec_wp$model <- 'Adjusted/Unmatched'


# SEX ATT
ame_3ni_unmatched_sex <- comparisons(model_unmatched3ni, variables="race_of_victim_black",
                                     by="sex_of_victim", newdata = subset(nibrs, race_of_victim_black == 1))|>
  summary()

ame_3ni_unmatched_sex$dataset <- 'NIBRS (1991-2020)'
ame_3ni_unmatched_sex$model <- 'Adjusted/Unmatched'

ame_3ni_unmatched_sex_pw <- comparisons(model_unmatched3ni, variables="race_of_victim_black",
                                        by="sex_of_victim", 
                                        hypothesis="pairwise", newdata = subset(nibrs, race_of_victim_black == 1))|>
  summary()

ame_3ni_unmatched_sex_pw$dataset <- 'NIBRS (1991-2020)'
ame_3ni_unmatched_sex_pw$model <- 'Adjusted/Unmatched'




###########  MAP, AFTER 1990 (SAME TIME WINDOW AS NIBRS)
map91<- filter(map, year2 > 1990) # create dataset keeping only murders occurred after 91 (included)

############## MAP (agency info + states)
m.out3r <- matchit(vic_race_black ~ vic_count + off_count + vic_sex+
                     age_cat+
                     decade+
                     agentype+
                     monthly_agency_overlap+
                     state,                        
                   data=map91,
                   method = "exact",
                   estimand="ATT")




bal.tab(m.out3r)

new.names2 <- c(state_Idaho = "State: Idaho",
                state_Utah = "State: Utah",
                agentype_city = "Agency: City",
                state_Montana = "State: Montana",
                agentype_state_police = "Agency: State Police",
                monthly_agency_overlap = "Monthly Agency Homicide Overlap",
                sex_of_victim_male = "Victim's Sex: Male",
                sex_of_victim_female = "Victim's Sex: Female",
                "state_New Hampshire" = "State: New Hampshire",
                state_Oregon ="State: Oregon",
                state_Colorado="State: Colorado",
                "state_South Dakota"="State: South Dakota",
                "state_West Virginia"="State: West Virginia",
                "state_North Dakota"="State: North Dakota",
                state_Washington="State: Washington",
                state_Vermont="State: Vermont",
                state_Arizona="State: Arizona",
                "state_New Mexico"="State: New Mexico",
                state_Hawaii="State: Hawaii",
                state_Oklahoma ="State: Oklahoma",
                "age_cat_21-25"="Victim's Age: 21-25",
                "age_cat_71-75"="Victim's Age: 71-75",
                "age_cat_76-80"="Victim's Age: 76-80",
                "age_cat_81-85"="Victim's Age: 81-85",
                state_Massachusetts="State: Massachusetts",
                "age_cat_61-65"="Victim's Age: 61-65",
                state_Michigan="State: Michigan",
                "age_cat_16-20"="Victim's Age: 16-20",
                "age_cat_66-70"="Victim's Age: 66-70",
                state_Iowa="State: Iowa",
                state_Nebraska="State: Nebraska",
                state_Ohio="State: Ohio",
                "age_cat_56-60"="Victim's Age: 56-60",
                agentype_other_agencies="Agency: Other Agencies",
                state_Kansas="State: Kansas",
                "state_South Carolina"="State: South Carolina",
                "age_cat_51-55"="Victim's Age: 51-55",
                "age_cat_26-30"="Victim's Age: 26-30",
                state_Maine="State: Maine",
                state_Kentucky="State: Kentucky",
                "age_cat_0-5"="Victim's Age: 0-5",
                total_victims="N of victims",
                "age_cat_46-50"="Victim's Age: 46-50",
                "age_cat_86-90"="Victim's Age: 86-90",
                agentype_tribal="Agency: Tribal",
                state_Virginia="State: Virginia",
                state_Texas="State: Texas",
                "age_cat_41-45"="Victim's Age: 41-45",
                total_offender_segments="N of offenders",
                state_Tennessee="State: Tennessee",
                "age_cat_91-95"="Victim's Age: 91-95",
                state_Wisconsin="State: Wisconsin",
                state_Missouri="State: Missouri",
                state_Georgia="State: Georgia",
                "state_North Carolina"="State: North Carolina",
                "state_Delaware"="State: Delaware",
                "age_cat_6-10"="Victim's Age: 6-10",
                state_Alabama="State: Alabama",
                "age_cat_31-35"="Victim's Age: 31-35",
                decade_90s="Decade: 1990s",
                "state_Rhode Island"="State: Rhode Island",
                state_Louisiana="State: Louisiana",
                decade_10s="Decade: 2010s",
                state_Nevada="State: Nevada",
                "age_cat_36-40"="Victim's Age: 36-40",
                state_Arkansas="State: Arkansas",
                state_Mississippi="State: Mississippi",
                state_Maryland="State: Maryland",
                state_Connecticut="State: Connecticut ",
                "age_cat_96-99"="Victim's Age: 96-99",
                state_Illinois="State: Illinois",
                decade_00s="Decade: 2000s",
                state_Wyoming="State: Wyoming",
                "state_New York"="State: New York",
                state_Minnesota="State: Minnesota",
                state_Indiana="State: Indiana",
                state_Pennsylvania="State: Pennsylvania",
                "agentype_covered-by another agency"="Agency: covered by oth. agency",
                "state_District Of Columbia"="State: DC ",
                "agentype_university or college"="Agency: University/college",
                "age_cat_11-15"="Victim's Age: 11-15",
                state_Florida="State: Florida",
                state_California="State: California",
                "agentype_Regional police"="Agency: Regional Police",
                "agentype_County police"="Agency: County Police",
                "agentype_Special police"="Agency: Special Police",
                "state_New Jersey"="State: New Jersey",
                "vic_sex_Unknown"="Victim's Sex: Unknown",
                vic_sex_Male="Victim's Sex: Male",
                vic_sex_Female="Victim's Sex: Female",
                agentype_Sheriff="Agency: Sheriff",
                "agentype_Municipal police"="Agency: Municipal Police",
                "agentype_Primary state LE"="Agency: Primary State LE",
                state_Alaska="State: Alaska",
                vic_count="N of victims",
                off_count="N of offenders"
)

l <- love.plot(m.out3r, stats = c("mean.diffs"),
          thresholds = c(m = .1, v = 2), abs = TRUE, 
          binary = "std",
          var.order = "unadjusted",
          var.names=new.names2,
          colors=c("sienna3", "deepskyblue3"),
          sample.names=c("Unmatched", "Matched"),
          title=NULL)

l + scale_x_continuous(labels = function(x) {
  ifelse(x < 1, gsub("^0", "", format(x, nsmall = 2)), format(x, nsmall = 2))
})

matched_dataset3r <- match.data(m.out3r)


# modeling adjusted and unmatched



# adjusted
model_matched3r <- glm(solved_yes ~ vic_race_black + 
                         vic_count + off_count + vic_sex+
                         age_cat + 
                         decade+
                         weapon+
                         circumstance+
                         agentype+	 
                         factor(monthly_agency_overlap)+
                         state, 
                       family=quasibinomial(), data = matched_dataset3r, weights=weights)


tidy(model_matched3r)

summ(model_matched3r, exp = TRUE, robust=TRUE, digits=3)


# GENERAL ATT
ame_3_91<-comparisons(model_matched3r, variables = "vic_race_black", 
                      vcov="HC3", wts="weights",
                      newdata = subset(matched_dataset3r, vic_race_black == 1)) |>
  summary()

ame_3_91$dataset <- 'MAP (1991-2020)'
ame_3_91$model <- 'Adjusted/Matched'


# simulation-based
sim_model_matched3r <- clarify::sim(model_matched3r)
est_sim_model_matched3r <- clarify::sim_ame(sim_model_matched3r, var="vic_race_black", subset=vic_race_black==1, contrast="diff", verbose=TRUE)
plot_est_sim_model_matched3r <- plot(est_sim_model_matched3r, color='red')+ylab("Densità")
summary(est_sim_model_matched3r)



# DECADE ATT
ame_3_91_dec<-comparisons(model_matched3r, variables = "vic_race_black", by="decade", 
                          vcov="HC3",
                          wts="weights", newdata = subset(matched_dataset3r, vic_race_black == 1)) |>
  summary()

ame_3_91_dec$dataset <- 'MAP (1991-2020)'
ame_3_91_dec$model <- 'Adjusted/Matched'


ame_3_91_dec_pw<-comparisons(model_matched3r, variables = "vic_race_black", by="decade", 
                             vcov="HC3",
                             wts="weights", newdata = subset(matched_dataset3r, vic_race_black == 1),
                             hypothesis="pairwise") |>
  summary()

ame_3_91_dec_pw$dataset <- 'MAP (1991-2020)'
ame_3_91_dec_pw$model <- 'Adjusted/Matched'



# SEX ATT
ame_3_91_sex<-comparisons(model_matched3r, variables = "vic_race_black", by="vic_sex", 
                          vcov="HC3",
                          wts="weights",
                          newdata = subset(matched_dataset3r, vic_race_black == 1)) |>
  summary()

ame_3_91_sex$dataset <- 'MAP (1991-2020)'
ame_3_91_sex$model <- 'Adjusted/Matched'



ame_3_91_sex_pw<-comparisons(model_matched3r, variables = "vic_race_black", by="vic_sex", 
                             vcov="HC3",
                             wts="weights",
                             newdata = subset(matched_dataset3r, vic_race_black == 1),
                             hypothesis="pairwise") |>
  summary()

ame_3_91_sex_pw$dataset <- 'MAP (1991-2020)'
ame_3_91_sex_pw$model <- 'Adjusted/Matched'



############ non matched
model_unmatched_3r <- glm(solved_yes ~ vic_race_black + 
                            vic_count + off_count + vic_sex+
                            age_cat + 
                            decade+
                            weapon+
                            circumstance+
                            agentype+	 
                            factor(monthly_agency_overlap)+
                            state, 
                          family=binomial, data = map91)

summ(model_unmatched_3r, exp = TRUE, robust=TRUE, digits=3)

# GENERAL ATT

ame_3_91_unmatched<-comparisons(model_unmatched_3r, variables = "vic_race_black", newdata = subset(map91, vic_race_black == 1)) |>
  summary()

ame_3_91_unmatched$dataset <- 'MAP (1991-2020)'
ame_3_91_unmatched$model <- 'Adjusted/Unmatched'

# simulation-based
sim_model_unmatched3r <- clarify::sim(model_unmatched_3r)
est_sim_model_unmatched3r <- clarify::sim_ame(sim_model_unmatched3r, var="vic_race_black", subset=vic_race_black==1, contrast="diff", verbose=TRUE)

plot_est_sim_model_unmatched3r <- plot(est_sim_model_unmatched3r, alpha=.8, fill="#FF6666")+ylab("Densità")+
  theme(strip.background = element_rect(fill = "yellow")) + theme_bw()
summary(est_sim_model_unmatched3r)




# DECADE ATT
ame_3_91_unmatched_dec<-comparisons(model_unmatched_3r, variables = "vic_race_black", by="decade", newdata = subset(map91, vic_race_black == 1)) |>
  summary()

ame_3_91_unmatched_dec$dataset <- 'MAP (1991-2020)'
ame_3_91_unmatched_dec$model <- 'Adjusted/Unmatched'



ame_3_91_unmatched_dec_pw<-comparisons(model_unmatched_3r, variables = "vic_race_black", 
                                       by="decade", newdata = subset(map91, vic_race_black == 1),
                                       hypothesis="pairwise") |>
  summary()

ame_3_91_unmatched_dec_pw$dataset <- 'MAP (1991-2020)'
ame_3_91_unmatched_dec_pw$model <- 'Adjusted/Unmatched'


# SEX ATT
ame_3_91_unmatched_sex<-comparisons(model_unmatched_3r, variables = "vic_race_black", by="vic_sex", newdata = subset(map91, vic_race_black == 1)) |>
  summary()

ame_3_91_unmatched_sex$dataset <- 'MAP (1991-2020)'
ame_3_91_unmatched_sex$model <- 'Adjusted/Unmatched'


ame_3_91_unmatched_sex_pw<-comparisons(model_unmatched_3r, variables = "vic_race_black", by="vic_sex", newdata = subset(map91, vic_race_black == 1),
                                       hypothesis="pairwise") |>
  summary()

ame_3_91_unmatched_sex_pw$dataset <- 'MAP (1991-2020)'
ame_3_91_unmatched_sex_pw$model <- 'Adjusted/Unmatched'




############### PLOTTING
# plotting (Figure on national AMEs)

ame_list <- do.call("rbind", list(ame_3_91, ame_3_91_unmatched,
                                  ame_3ni, ame_3ni_unmatched ))
ame_list <- ame_list[order(ame_list$dataset),]
ame_list$dataset<- fct_inorder(ame_list$dataset)
ame_list$model = factor(ame_list$model, levels=c("Adjusted/Matched", "Adjusted/Unmatched"))
# plot
dodge <- position_dodge(width=0.5)
ame_plot <- ggplot(ame_list, x=factor(ame_list$dataset, levels=unique(ame_list$dataset)), y=estimate, color=model, aes(ymin=-0.06,ymax=0.06))+
  geom_point(aes(x=factor(ame_list$dataset, levels=unique(ame_list$dataset)), y=estimate,shape=model, color=model), size=2.5, position=dodge)+
  scale_y_continuous(breaks =c("-.06"=-0.06, "-.04"=-0.04,"-.02"=-0.02, "0"=0,
                               0, ".02"=0.02, ".04"=0.04, ".06"=0.06))+
  geom_errorbar(aes(x=dataset, ymin=conf.low, ymax=conf.high, group=model, color=model), width=0.1,position=dodge, name="Model", labels = c("Adjusted/Matched", "Adjusted/Unmatched"))+
  scale_color_manual(values=c("Adjusted/Matched"="red","Adjusted/Unmatched"="blue"), name="Model", labels = c("Adjusted/Matched", "Adjusted/Unmatched"))+
  scale_shape_manual(values=c("Adjusted/Matched"=19,"Adjusted/Unmatched"= 15), name="Model", labels = c("Adjusted/Matched", "Adjusted/Unmatched" ))+
  labs(x= "Dataset", y = "Average Marginal Effect (95% C.I.)")+
  theme_bw()+
  geom_hline(aes(yintercept=0), colour="#990000", linetype="dashed")+
  coord_flip()+theme(axis.text.y = element_text(size=12)) + 
  theme(axis.text.x = element_text(size=12), axis.title=element_text(size=16))+
  #theme(text = element_text(family = "Times New Roman"))+
  theme(legend.title=element_text(size=14), legend.key.size=unit(0.5, 'cm'), legend.text = element_text(size=14))+theme(legend.position="bottom")+
  guides(colour=guide_legend(nrow=1,byrow=TRUE))

# customized plot NIBRS matched
df_sim_model_matched3ni <- as.data.frame(est_sim_model_matched3ni)
df_sim_model_matched3ni$ID <- seq.int(nrow(df_sim_model_matched3ni))
df_sim_model_matched3ni_long <- melt(setDT(df_sim_model_matched3ni), id.vars="ID", variable.name="quantity")


plot_sim_model_matched3ni <- ggplot(aes(x=value), data=df_sim_model_matched3ni_long)

plot_sim_model_matched3ni<-plot_sim_model_matched3ni+geom_density(aes(fill=quantity), alpha=0.4)+
  scale_fill_manual( values = c("orange2","slategrey", "darkturquoise"))+
  geom_hline(yintercept=0)+
  facet_wrap(. ~ quantity, scales='free', labeller=as_labeller(c(`E[Y(0)]`="E[Clearance|Victim: non-Black]", 
                                                                 `E[Y(1)]` = "E[Clearance|Victim: Black]", 
                                                                 `Diff`="Difference (AME)")))+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="E[Y(0)]"), aes(xintercept=0.5231), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="E[Y(0)]"), aes(xintercept=0.5343), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="E[Y(0)]"), aes(xintercept=0.5289), colour="black", linetype="solid", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="E[Y(1)]"), aes(xintercept=0.4904), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="E[Y(1)]"), aes(xintercept=0.4990), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="E[Y(1)]"), aes(xintercept=0.4946), colour="black", linetype="solid", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="Diff"), aes(xintercept=-0.0412), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="Diff"), aes(xintercept=-0.0270), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="Diff"), aes(xintercept=-0.0343), colour="black", linetype="solid", linewidth=1)+
  ylab('Density')+xlab("Estimates (NIBRS Adjusted/Matched)")+
  theme_bw()+theme(legend.position = "none", strip.background = element_rect(colour="black",
                                                                             fill="white"))+
  scale_x_continuous(labels = function(x) {
    ifelse(x == 0, "0", 
           ifelse(x > 0, gsub("^0", "", format(x, nsmall = 3)), 
                  gsub("^-0", "-.", gsub("^-0\\.", "-.", format(x, nsmall = 3)))
           )
    )
  })

# customized plot NIBRS unmatched
df_sim_model_unmatched3ni <- as.data.frame(est_sim_model_unmatched3ni)
df_sim_model_unmatched3ni$ID <- seq.int(nrow(df_sim_model_unmatched3ni))
df_sim_model_unmatched3ni_long <- melt(setDT(df_sim_model_unmatched3ni), id.vars="ID", variable.name="quantity")


plot_sim_model_unmatched3ni <- ggplot(aes(x=value), data=df_sim_model_unmatched3ni_long)

plot_sim_model_unmatched3ni<-plot_sim_model_unmatched3ni+geom_density(aes(fill=quantity), alpha=0.4)+
  scale_fill_manual( values = c("orange2","slategrey", "darkturquoise"))+
  geom_hline(yintercept=0)+
  facet_wrap(. ~ quantity, scales='free', labeller=as_labeller(c(`E[Y(0)]`="E[Clearance|Victim: non-Black]", 
                                                                 `E[Y(1)]` = "E[Clearance|Victim: Black]", 
                                                                 `Diff`="Difference (AME)")))+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="E[Y(0)]"), aes(xintercept=0.5578), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="E[Y(0)]"), aes(xintercept=0.5699), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="E[Y(0)]"), aes(xintercept=0.56419), colour="black", linetype="solid", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="E[Y(1)]"), aes(xintercept=0.5188), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="E[Y(1)]"), aes(xintercept=0.5268), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="E[Y(1)]"), aes(xintercept=0.5228), colour="black", linetype="solid", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="Diff"), aes(xintercept=-0.0482), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="Diff"), aes(xintercept=-0.0337), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3ni_long, quantity=="Diff"), aes(xintercept=-0.0413), colour="black", linetype="solid", linewidth=1)+
  ylab('Density')+xlab("Estimates (NIBRS Adjusted/Unmatched)")+
  theme_bw()+theme(legend.position = "none", strip.background = element_rect(colour="black",
                                                                             fill="white"))+
  scale_x_continuous(labels = function(x) {
    ifelse(x == 0, "0", 
           ifelse(x > 0, gsub("^0", "", format(x, nsmall = 3)), 
                  gsub("^-0", "-.", gsub("^-0\\.", "-.", format(x, nsmall = 3)))
           )
    )
  })

# customised plot MAP matched
df_sim_model_matched3r <- as.data.frame(est_sim_model_matched3r)
df_sim_model_matched3r$ID <- seq.int(nrow(df_sim_model_matched3r))
df_sim_model_matched3r_long <- melt(setDT(df_sim_model_matched3r), id.vars="ID", variable.name="quantity")


plot_sim_model_matched3r <- ggplot(aes(x=value), data=df_sim_model_matched3r_long)

plot_sim_model_matched3r<-plot_sim_model_matched3r+geom_density(aes(fill=quantity), alpha=0.4)+
  scale_fill_manual( values = c("orange2","slategrey", "darkturquoise"))+
  geom_hline(yintercept=0)+
  facet_wrap(. ~ quantity, labeller=as_labeller(c(`E[Y(0)]`="E[Clearance|Victim: non-Black]", 
                                                  `E[Y(1)]` = "E[Clearance|Victim: Black]", 
                                                  `Diff`="Difference (AME)")), scales='free')+
  geom_vline(data=subset(df_sim_model_matched3r_long, quantity=="E[Y(0)]"), aes(xintercept=0.6426), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3r_long, quantity=="E[Y(0)]"), aes(xintercept=0.6661), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3r_long, quantity=="E[Y(0)]"), aes(xintercept=0.6555), colour="black", linetype="solid", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3r_long, quantity=="E[Y(1)]"), aes(xintercept=0.5989), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3r_long, quantity=="E[Y(1)]"), aes(xintercept=0.6205), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3r_long, quantity=="E[Y(1)]"), aes(xintercept=0.6107), colour="black", linetype="solid", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3r_long, quantity=="Diff"), aes(xintercept=-0.0598), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3r_long, quantity=="Diff"), aes(xintercept=-0.0294), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_matched3r_long, quantity=="Diff"), aes(xintercept=-0.0448), colour="black", linetype="solid", linewidth=1)+
  ylab('Density')+xlab("Estimates (MAP Adjusted/Matched)")+
  theme_bw()+theme(legend.position = "none", strip.background = element_rect(colour="black",
                                                                             fill="white"))+
  scale_x_continuous(labels = function(x) {
    ifelse(x == 0, "0", 
           ifelse(x > 0, gsub("^0", "", format(x, nsmall = 3)), 
                  gsub("^-0", "-.", gsub("^-0\\.", "-.", format(x, nsmall = 3)))
           )
    )
  })


# customised plot MAP unmatched
df_sim_model_unmatched3r <- as.data.frame(est_sim_model_unmatched3r)
df_sim_model_unmatched3r$ID <- seq.int(nrow(df_sim_model_unmatched3r))
df_sim_model_unmatched3r_long <- melt(setDT(df_sim_model_unmatched3r), id.vars="ID", variable.name="quantity")


plot_sim_model_unmatched3r <- ggplot(aes(x=value), data=df_sim_model_unmatched3r_long)

plot_sim_model_unmatched3r<-plot_sim_model_unmatched3r+geom_density(aes(fill=quantity), alpha=0.4)+
  scale_fill_manual( values = c("orange2","slategrey", "darkturquoise"))+
  geom_hline(yintercept=0)+
  facet_wrap(. ~ quantity, labeller=as_labeller(c(`E[Y(0)]`="E[Clearance|Victim: non-Black]", 
                                                  `E[Y(1)]` = "E[Clearance|Victim: Black]", 
                                                  `Diff`="Difference (AME)")), scales='free')+
  geom_vline(data=subset(df_sim_model_unmatched3r_long, quantity=="E[Y(0)]"), aes(xintercept=0.6741), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_unmatched3r_long, quantity=="E[Y(0)]"), aes(xintercept=0.6700), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_unmatched3r_long, quantity=="E[Y(0)]"), aes(xintercept=0.6720), colour="black", linetype="solid", linewidth=1)+
  geom_vline(data=subset(df_sim_model_unmatched3r_long, quantity=="E[Y(1)]"), aes(xintercept=0.6247), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_unmatched3r_long, quantity=="E[Y(1)]"), aes(xintercept=0.6217), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_unmatched3r_long, quantity=="E[Y(1)]"), aes(xintercept=0.6232), colour="black", linetype="solid", linewidth=1)+
  geom_vline(data=subset(df_sim_model_unmatched3r_long, quantity=="Diff"), aes(xintercept=-0.0512), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_unmatched3r_long, quantity=="Diff"), aes(xintercept=-0.0462), colour="pink", linetype="dashed", linewidth=1)+
  geom_vline(data=subset(df_sim_model_unmatched3r_long, quantity=="Diff"), aes(xintercept=-0.0488), colour="black", linetype="solid", linewidth=1)+
  ylab('Density')+xlab("Estimates (MAP Adjusted/Unmatched)")+
  theme_bw()+theme(legend.position = "none", strip.background = element_rect(colour="black",
                                                                             fill="white"))+
  scale_x_continuous(labels = function(x) {
    ifelse(x == 0, "0", 
           ifelse(x > 0, gsub("^0", "", format(x, nsmall = 3)), 
                  gsub("^-0", "-.", gsub("^-0\\.", "-.", format(x, nsmall = 3)))
           )
    )
  })


ggarrange(ame_plot, 
          ggarrange(plot_sim_model_matched3r,
                    plot_sim_model_unmatched3r,
                    plot_sim_model_matched3ni,
                    plot_sim_model_unmatched3ni), labels=c("A", "B"), heights=c(0.8,1),
          nrow=2, ncol=1)

ggarrange(ame_plot, 
          ggarrange(plot_sim_model_matched3r,
                    plot_sim_model_unmatched3r),
          ggarrange(plot_sim_model_matched3ni,
                    plot_sim_model_unmatched3ni),
          labels=c("A", "B", ""), heights=c(0.8,1,1),
          nrow=3, ncol=1)

ggarrange(ame_plot, plot_sim_model_matched3r,plot_sim_model_unmatched3r, plot_sim_model_matched3ni, plot_sim_model_unmatched3ni,
          labels=c("A", "B", "", "", ""), heights=c(1.3,0.8,0.8,0.8,0.8),
          nrow=5, ncol=1)

####################### sex plot

### descriptive

## MAP (Full)

vicsex_map <- map91 %>% group_by(vic_sex, vic_race_black) %>% summarise(Count = n(),
                                                                        Percentage=n()/nrow(.), 
                                                                        Mean = mean(solved_yes))

vicsex_map_w <- reshape(as.data.frame(vicsex_map), idvar="vic_sex",
                        v.names=c("Count", "Percentage", "Mean"),
                        timevar="vic_race_black",
                        direction="wide")

xtable(vicsex_map)

vicsex_map$vic_race_black <- as.factor(vicsex_map$vic_race_black)


vicsex_map_long <- vicsex_map %>%  pivot_longer(cols=c('Mean', 'Percentage'), names_to="Type", values_to='ratios')
vicsex_map_long$dataset <- 'MAP (Full)'
vicsex_map_plot <- ggplot(vicsex_map_long)+
  geom_point(vicsex_map_long, mapping=aes(y=vic_sex, x=ratios, shape=Type, color=vic_race_black),size=4, alpha=0.8)+
  theme_bw()+ylab("Victim's Sex")+xlab("Ratio")+scale_color_manual(values=c("royalblue1", "indianred3"), labels=c("Not Black", "Black"))+
  scale_shape_manual(values=c(4,3), labels=c("Ratio of Solved", "Ratio of Cases"))+labs(shape="Type of Ratio", colour="Victim's Race")+
  xlim(0,1)+theme(axis.text=element_text(size=12),
                  axis.title=element_text(size=16))

## MAP (Matched)

vicsex_map_m <- matched_dataset3r %>% group_by(vic_sex, vic_race_black) %>% summarise(Count = n(),
                                                                                      Percentage=n()/nrow(.), 
                                                                                      Mean = mean(solved_yes))



vicsex_map_m$vic_race_black <- as.factor(vicsex_map_m$vic_race_black)


vicsex_map_long_m <- vicsex_map_m %>%  pivot_longer(cols=c('Mean', 'Percentage'), names_to="Type", values_to='ratios')
vicsex_map_long_m$dataset <- 'MAP (Matched)'
vicsex_map_plot_m <- ggplot(vicsex_map_long_m)+
  geom_point(vicsex_map_long_m, mapping=aes(y=vic_sex, x=ratios, shape=Type, color=vic_race_black),size=4, alpha=0.8)+
  theme_bw()+ylab("Victim's Sex")+xlab("Ratio")+scale_color_manual(values=c("royalblue1", "indianred3"), labels=c("Not Black", "Black"))+
  scale_shape_manual(values=c(4,3), labels=c("Ratio of Solved", "Ratio of Cases"))+labs(shape="Type of Ratio", colour="Victim's Race")+
  xlim(0,1)+theme(axis.text=element_text(size=12),
                  axis.title=element_text(size=16))

## NIBRS (Full)
vicsex_nibrs <- nibrs %>% group_by(sex_of_victim, race_of_victim_black) %>% summarise(Count = n(),
                                                                                      Percentage=n()/nrow(.), 
                                                                                      Mean = mean(solved))



vicsex_nibrs_long <- vicsex_nibrs %>% pivot_longer(cols=c('Mean', 'Percentage'), names_to="Type", values_to="ratios")
vicsex_nibrs_long$dataset <- 'NIBRS (Full)'

vicsex_nibrs_long <- rename(vicsex_nibrs_long, vic_race_black = race_of_victim_black)
vicsex_nibrs_long <- rename(vicsex_nibrs_long, vic_sex = sex_of_victim)
vicsex_nibrs_long$vic_sex[vicsex_nibrs_long$vic_sex=='female'] <- 'Female'
vicsex_nibrs_long$vic_sex[vicsex_nibrs_long$vic_sex=='male'] <- 'Male'
vicsex_nibrs_long$vic_sex[vicsex_nibrs_long$vic_sex=='unknown'] <- 'Unknown'

## NIBRS (Matched)

## NIBRS (Full)
vicsex_nibrs_m <- matched_dataset3ni %>% group_by(sex_of_victim, race_of_victim_black) %>% summarise(Count = n(),
                                                                                                     Percentage=n()/nrow(.), 
                                                                                                     Mean = mean(solved))



vicsex_nibrs_long_m <- vicsex_nibrs_m %>% pivot_longer(cols=c('Mean', 'Percentage'), names_to="Type", values_to="ratios")
vicsex_nibrs_long_m$dataset <- 'NIBRS (Matched)'

vicsex_nibrs_long_m <- rename(vicsex_nibrs_long_m, vic_race_black = race_of_victim_black)
vicsex_nibrs_long_m <- rename(vicsex_nibrs_long_m, vic_sex = sex_of_victim)
vicsex_nibrs_long_m$vic_sex[vicsex_nibrs_long_m$vic_sex=='female'] <- 'Female'
vicsex_nibrs_long_m$vic_sex[vicsex_nibrs_long_m$vic_sex=='male'] <- 'Male'
vicsex_nibrs_long_m$vic_sex[vicsex_nibrs_long_m$vic_sex=='unknown'] <- 'Unknown'

# combine map + nibrs

vicsex_long <- do.call("rbind", list(vicsex_map_long, vicsex_map_long_m, vicsex_nibrs_long, vicsex_nibrs_long_m))

vicsex_plot <- ggplot(vicsex_long)+
  geom_point(vicsex_long, mapping=aes(y=vic_sex, x=ratios, shape=Type, color=as.factor(vic_race_black)),size=4, alpha=1)+
  theme_bw()+ylab("Victim's Sex")+xlab("Ratio")+scale_color_manual(values=c("orange2", "slategrey"), labels=c("Not Black", "Black"))+
  scale_shape_manual(values=c(4,3), labels=c("Ratio of Solved", "Ratio of Cases"))+labs(shape="Type of Ratio", colour="Victim's Race")+
  scale_y_discrete(labels=c("Female", "Male", "Unknown"))+
  xlim(0,1)+theme(axis.text=element_text(size=12),
                  axis.title=element_text(size=16))+facet_wrap(. ~ dataset)+theme(legend.position = "bottom", strip.background = element_rect(colour="black",
                                                                                                                                              fill="white"))+
  scale_x_continuous(limits = c(-0.05, 1.05), expand = c(0, 0),
                     labels = function(x) {
                       ifelse(x == 0, "0", 
                              ifelse(x > 0, gsub("^0", "", format(x, nsmall = 2)), 
                                     gsub("^-0", "-.", gsub("^-0\\.", "-.", format(x, nsmall = 2)))
                              )
                       )
                     })



vicsex_nibrs_plot <- ggplot(vicsex_nibrs_long)+
  geom_point(vicsex_nibrs_long, mapping=aes(y=sex_of_victim, x=ratios, shape=Type, color=as.factor(race_of_victim_black)),size=4, alpha=1)+
  theme_bw()+ylab("Victim's Sex")+xlab("Ratio")+scale_color_manual(values=c("royalblue1", "indianred3"), labels=c("Not Black", "Black"))+
  scale_shape_manual(values=c(4,3), labels=c("Ratio of Solved", "Ratio of Cases"))+labs(shape="Type of Ratio", colour="Victim's Race")+
  scale_y_discrete(labels=c("Female", "Male", "Unknown"))+
  xlim(0,1)+theme(axis.text=element_text(size=12),
                  axis.title=element_text(size=16))





## group-ame
ame_3ni_sex <- rename(ame_3ni_sex, c('vic_sex'='sex_of_victim'))
ame_3ni_unmatched_sex <- rename(ame_3ni_unmatched_sex, c('vic_sex'='sex_of_victim'))

ame_sex_list <- do.call("rbind", list(ame_3_91_sex, ame_3_91_unmatched_sex,
                                      ame_3ni_sex, ame_3ni_unmatched_sex))

ame_sex_list <- ame_sex_list[order(ame_sex_list$dataset),]
ame_sex_list$dataset <- fct_inorder(ame_sex_list$dataset)
ame_sex_list$model = factor(ame_sex_list$model, levels=c("Adjusted/Matched", "Adjusted/Unmatched"))
#capitalize sex for consistency and easier plotting
ame_sex_list$vic_sex <- str_to_title(ame_sex_list$vic_sex)


dodge <- position_dodge(width=0.4)
ame_sex_plot <- ggplot(ame_sex_list, x=factor(ame_sex_list$vic_sex, levels=unique(ame_sex_list$vic_sex)), y=estimate, color=model, aes(ymin=-0.08,ymax=0.08))+
  geom_point(aes(x=factor(ame_sex_list$vic_sex, levels=unique(ame_sex_list$vic_sex)), y=estimate,shape=model, color=model), size=2.5, position=dodge)+
  scale_y_continuous(breaks =c("-.10"=-0.10, "-.08"=-0.08,"-.06"=-0.06, "-.04"=-0.04,"-.02"=-0.02, "0"=0, 
                               ".02"=0.02, ".04"=0.04, ".06"=0.06, ".08"=0.08, ".10"=0.10))+
  geom_errorbar(aes(x=vic_sex, ymin=conf.low, ymax=conf.high, group=model, color=model), width=0.1,position=dodge, name="Model", labels = c("Adjusted/Matched", "Adjusted/Unmatched"))+
  scale_color_manual(values=c("Adjusted/Matched"="red","Adjusted/Unmatched"="blue"), name="Model", labels = c("Adjusted/Matched", "Adjusted/Unmatched"))+
  scale_shape_manual(values=c("Adjusted/Matched"=19,"Adjusted/Unmatched"= 15), name="Model", labels = c("Adjusted/Matched", "Adjusted/Unmatched" ))+
  labs(x= "Victim's Sex", y = "Average Marginal Effect (95% C.I.)")+
  theme_bw()+
  geom_hline(aes(yintercept=0), colour="#990000", linetype="dashed")+
  coord_flip()+theme(axis.text.y = element_text(size=12)) + 
  theme(axis.text.x = element_text(size=12), axis.title=element_text(size=16))+
  #theme(text = element_text(family = "Times New Roman"))+
  theme(legend.title=element_text(size=13), legend.key.size=unit(0.1, 'cm'), legend.text = element_text(size=10))+theme(legend.position="bottom")+
  guides(colour=guide_legend(nrow=1,byrow=TRUE))

ame_sex_plot <- ame_sex_plot + facet_grid(.~ factor(dataset), labeller=as_labeller(c(`MAP (1991-2020)`="MAP", 
                                                                                     `NIBRS (1991-2020)` = "NIBRS")))+ theme(strip.background = element_rect(fill = "white"),
                                                                                                                             strip.text = element_text(color = "black"))

# combine
ggarrange(vicsex_plot, ame_sex_plot, ncol=1, labels=c("A", "B"))

####################### decade plot

### descriptive


# MAP (Full)

decade_map <- map91 %>% group_by(decade, vic_race_black) %>% summarise(Count = n(),
                                                                       Percentage=n()/nrow(.), 
                                                                       Mean = mean(solved_yes))

xtable(decade_map)

decade_map_long <- decade_map %>%  pivot_longer(cols=c('Mean', 'Percentage'), names_to="Type", values_to='ratios')

decade_map_long$decade <- factor(decade_map_long$decade, levels=c("90s", "00s", "10s"))
decade_map_long$dataset <- 'MAP (Full)'


# MAP (Matched)
decade_map_m <- matched_dataset3r %>% group_by(decade, vic_race_black) %>% summarise(Count = n(),
                                                                                     Percentage=n()/nrow(.), 
                                                                                     Mean = mean(solved_yes))


decade_map_long_m <- decade_map_m %>%  pivot_longer(cols=c('Mean', 'Percentage'), names_to="Type", values_to='ratios')

decade_map_long_m$decade <- factor(decade_map_long_m$decade, levels=c("90s", "00s", "10s"))
decade_map_long_m$dataset <- 'MAP (Matched)'

## NIBRS (Full)
decade_nibrs <- nibrs %>% group_by(decade, race_of_victim_black) %>% summarise(Count = n(),
                                                                               Percentage=n()/nrow(.), 
                                                                               Mean = mean(solved))

xtable(decade_nibrs)

decade_nibrs_long <- decade_nibrs %>% pivot_longer(cols=c('Mean', 'Percentage'), names_to="Type", values_to="ratios")

decade_nibrs_long$decade <- factor(decade_nibrs_long$decade, levels=c("90s", "00s", "10s"))
decade_nibrs_long$dataset<-'NIBRS (Full)'
decade_nibrs_long <- rename(decade_nibrs_long, vic_race_black = race_of_victim_black)

## NIBRS (Matched)
decade_nibrs_m <- matched_dataset3ni %>% group_by(decade, race_of_victim_black) %>% summarise(Count = n(),
                                                                                              Percentage=n()/nrow(.), 
                                                                                              Mean = mean(solved))

xtable(decade_nibrs)

decade_nibrs_long_m <- decade_nibrs_m %>% pivot_longer(cols=c('Mean', 'Percentage'), names_to="Type", values_to="ratios")

decade_nibrs_long_m$decade <- factor(decade_nibrs_long_m$decade, levels=c("90s", "00s", "10s"))
decade_nibrs_long_m$dataset<-'NIBRS (Matched)'
decade_nibrs_long_m <- rename(decade_nibrs_long_m, vic_race_black = race_of_victim_black)

# combine map + nibrs

decade_long <- do.call("rbind", list(decade_map_long, decade_map_long_m, decade_nibrs_long, decade_nibrs_long_m))

decade_plot <- ggplot(decade_long)+
  geom_point(decade_long, mapping=aes(y=decade, x=ratios, shape=Type, color=as.factor(vic_race_black)),size=4, alpha=1)+
  theme_bw()+ylab("Decade")+xlab("Ratio")+scale_color_manual(values=c("orange2", "slategrey"), labels=c("Not Black", "Black"))+
  scale_shape_manual(values=c(4,3), labels=c("Ratio of Solved", "Ratio of Cases"))+labs(shape="Type of Ratio", colour="Victim's Race")+
  scale_y_discrete(labels=c("1990s", "2000s", "2010s"))+
  xlim(0,1)+theme(axis.text=element_text(size=12),
                  axis.title=element_text(size=16))+facet_wrap(. ~ dataset)+theme(legend.position = "bottom", strip.background = element_rect(colour="black",
                                                                                                                                              fill="white"))+
  scale_x_continuous(limits = c(-0.05, 1.05), expand = c(0, 0),
                     labels = function(x) {
                       ifelse(x == 0, "0", 
                              ifelse(x > 0, gsub("^0", "", format(x, nsmall = 2)), 
                                     gsub("^-0", "-.", gsub("^-0\\.", "-.", format(x, nsmall = 2)))
                              )
                       )
                     })





### group-ame
ame_decade_list <- do.call("rbind", list(ame_3_91_dec, ame_3_91_unmatched_dec,
                                         ame_3ni_dec, ame_3ni_unmatched_dec))

ame_decade_list <- ame_decade_list[order(ame_decade_list$dataset),]
ame_decade_list$dataset <- fct_inorder(ame_decade_list$dataset)
ame_decade_list$model = factor(ame_decade_list$model, levels=c("Adjusted/Matched", "Adjusted/Unmatched"))

dodge <- position_dodge(width=0.4)


# Assuming your data frame is named ame_decade_list
ame_decade_list <- ame_decade_list %>%
  arrange(
    ifelse(decade == "90s", 1,
           ifelse(decade == "00s", 2,
                  ifelse(decade == "10s", 3, NA)))
  )

# Assuming your data frame is named ame_decade_list
ame_decade_list <- ame_decade_list %>%
  mutate(decade = ifelse(decade == "90s", "1990s",
                         ifelse(decade == "00s", "2000s",
                                ifelse(decade == "10s", "2010s", as.character(decade))))
  )

ame_decade_list$decade <- as.character(ame_decade_list$decade)


ame_decade_list$decade <- factor(ame_decade_list$decade, levels=c("1990s", "2000s", "2010s"))


ame_decade_plot <- ggplot(ame_decade_list, x=factor(ame_decade_list$decade, levels=unique(ame_decade_list$decade)), y=estimate, color=model, aes(ymin=-0.06,ymax=0.06))+
  geom_point(aes(x=factor(ame_decade_list$decade, levels=unique(ame_decade_list$decade)), y=estimate,shape=model, color=model), size=2.5, position=dodge)+
  scale_y_continuous(breaks =c("-.06"=-0.06, "-.04"=-0.04,"-.02"=-0.02, "0"=0, ".02"=0.02, ".04"=0.04, ".06"=0.06))+
  geom_errorbar(aes(x=decade, ymin=conf.low, ymax=conf.high, group=model, color=model), width=0.1,position=dodge, name="Model", labels = c("Adjusted/Matched", "Adjusted/Unmatched"))+
  scale_color_manual(values=c("Adjusted/Matched"="red","Adjusted/Unmatched"="blue"), name="Model", labels = c("Adjusted/Matched", "Adjusted/Unmatched"))+
  scale_shape_manual(values=c("Adjusted/Matched"=19,"Adjusted/Unmatched"= 15), name="Model", labels = c("Adjusted/Matched", "Adjusted/Unmatched" ))+
  labs(x= "Decade", y = "Average Marginal Effect (95% C.I.)")+
  theme_bw()+
  geom_hline(aes(yintercept=0), colour="#990000", linetype="dashed")+
  coord_flip()+theme(axis.text.y = element_text(size=12))+
  theme(axis.text.x = element_text(size=12), axis.title=element_text(size=16))+
  #theme(text = element_text(family = "Times New Roman"))+
  theme(legend.title=element_text(size=13), legend.key.size=unit(0.1, 'cm'), legend.text = element_text(size=10))+theme(legend.position="bottom")+
  guides(colour=guide_legend(nrow=1,byrow=TRUE))




ame_decade_plot <- ame_decade_plot + facet_grid(.~ factor(dataset), labeller=as_labeller(c(`MAP (1991-2020)`="MAP", 
                                                                                           `NIBRS (1991-2020)` = "NIBRS")))+theme(strip.background = element_rect(fill = "white"),
                                                                                                                                  strip.text = element_text(color = "black"))



# combine
ggarrange(decade_plot, ame_decade_plot, ncol=1, labels=c("A", "B"))










